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Abstract 

The Digital Orrery has been used to perform an in- 
tegration of the motion of the outer planets for 845 million 
years. This integration indicates that the long-term motion 
of the planet Pluto is chaotic. Nearby trajectories diverge 
exponentially with an e-folding time of only about 20 mil- 
lion years. 

The determination of the stability of the Solar System is one of the old- 
est problems in dynamical astronomy, but despite considerable attention all 
attempts to prove the stability of the system have failed. Arnold has shown 
that a large proportion of possible solar systems are stable if the masses, and 
orbital eccentricities and inclinations of the planets are sufficiently small [1]. 
The actual Solar System, however, does not meet the stringent requirements 
of the proof. Certainly, the great age of the Solar System suggests a high 
level of stability, but the nature of the long-term motion remains undeter- 
mined. The apparent analytical complexity of the problem of determining 
the stability of the solar system has led us to investigate the stability by 
means of numerical models. We have investigated the long-term stability of 
the Solar System through an 845 million year numerical integration of the 
outer planets with the Digital Orrery [2], a special purpose computer for 
studying planetary motion. 

Our earlier integrations [3] showed surprisingly long periods in the mo- 
tion of Pluto, of order 137 million years, and longer periods or a possible 



secular drift in the inclination of Pluto. The motion of Pluto is involved in a 
surprisingly large number of resonances. The libration of the resonant argu- 
ment permits the orbit of Pluto to cross the orbit of Neptune without having 
close encounters [4]. Williams and Benson [5] also found that the argument 
of perihelion was oscillating about x/2. Thus the perihelion of Pluto remains 
far from the line of intersection of the orbital planes of Pluto and Neptune. 
Williams and Benson noted that two other resonance conditions were pos- 
sibly satisfied by the motion of Pluto. Using the frequencies determined in 
our 200 million year integrations, we find that these two resonances are sat- 
isfied to the precision of our measurements. In addition, the 137 million year 
components that we observed in the elements of Pluto are associated with 
yet another commensurability. In Hamiltonian systems resonances are asso- 
ciated with both enhanced stability and instability. Without the libration 
of the resonant argument or the argument of perihelion a Neptune crossing 
orbit would be highly unstable. Nevertheless, resonances are also associated 
with most prominent zones of instability. This is particularly true if several 
resonances strongly affect the motion. 

The very long periods and the abundance of resonances in the motion 
of Pluto compelled us to pursue longer integrations. Our new integration 
indicates that the long-term motion of the planet Pluto, and by implication 
the motion of the whole Solar System, is weakly unstable. Exponential di- 
vergence of nearby trajectories indicates that the motion of Pluto is chaotic 
with the largest Lyapunov exponent estimated to be 10~ 7 ' 3 yr -1 . 

Our Numerical Experiment 

For many years the longest direct integration of the outer planets was the 
one million year integration of Cohen, Hubbard, and Oesterwinter [6]. Re- 
cently several new integrations of the outer planets have been performed [7,3,8]. 
The longest was our set of 200 million year integrations. Our new integra- 
tion is significantly longer and more accurate than all previously reported 
long-term integrations. 

In our new integration of the motion of the outer planets the masses 
and initial conditions were the same as those used in our 200 million year 
integrations of the outer planets. The planet Pluto is taken to be a zero- 
mass test particle. We continue to neglect the effects of the inner planets, 
the mass lost by the Sun due to electromagnetic radiation and solar wind, 



and general relativity. The most serious limitation of our integration is our 
ignorance of the true masses and initial conditions. We believe that our model 
is close enough to the actual solar sytem that its study sheds light on the 
question of stability of the solar system. To draw more rigorous conclusions 
the sensitivity of our conclusions to the uncertainties in masses and initial 
conditions, and to unmodeled effects must be determined. 

Our earlier integrations were limited to 100 million years forward and 
backward in time because of the accumulation of error, which was most 
seriously manifested in an accumulated longitude error of Jupiter of order 
50 degrees. In our new integrations we continue to use the 12 th order Stormer 
predictor, but a judicious choice of stepsize has reduced the numerical errors 
by several orders of magnitude. In all of our integrations the error in energy 
of the system varies nearly linearly with time. In the regime where neither 
roundoff nor truncation error is dominant the slope of the energy error as a 
function of time depends on stepsize in a complicated way. For some stepsizes 
the energy error has a positive slope; for others the slope is negative. This 
suggests that there might be special stepsizes for which there is no linear 
growth of energy error. By a series of numerical experiments we indeed 
found that there are values of the stepsize where the slope of the linear 
trend of energy vanishes. The special stepsizes become better defined as the 
integration interval of the experiments is increased. 

We chose our stepsize on the basis of a dozen 3 million year integrations, 
and numerous shorter integrations. For our new long integration we chose 
the stepsize to be 32.7 days. This seemingly innocuous change from a stepsize 
near 40 days dramatically reduces the slope of the energy error, by roughly 
three orders of magnitude. If the numerical integration were truncation-error 
dominated, for which the accumulated error is proportional to the h n , where 
h is the stepsize and n is the order of the integrator, then this reduction 
of stepsize would only improve the accumulated error by about a factor of 
about 10. 

In our new integration the relative energy error (energy minus initial 
energy divided by the magnitude of the initial energy) over 845 million years 
is —2.6 x 10" 10 ; the growth of the relative energy error is still very nearly linear 
with a slope of — 3.0 x 10~ 19 per year. By comparison the rate of growth of 
the relative energy error in our 200 million year integrations was 1.8 x 10" 16 
per year. The errors in other integrations of the outer Solar System were 



comparable to the errors in our 200 million year integrations. The rate of 
growth of energy error in the one million year integration of Cohen, Hubbard, 
and Oesterwinter was 2.4 x 10~ 16 per year. For the 6 million year integration 
of Kinoshita and Nakai the relative energy error was approximately 5 x 10~ 16 
per year. For the LONGSTOP integration the growth of relative energy (as 
defined in this paper) was -2.5 x 10" 16 per year. Thus the rate of growth 
of energy error in the integration reported in this paper is smaller than all 
previous long-term integrations of the outer planets by a factor of about 600. 

We verified that this improvement in energy conservation was reflected in 
a corresponding improvement in position and velocity errors by integrating 
the outer planets forward 3 million years and then backward to recover the 
initial conditions, over a range of stepsizes. For the chosen stepsize of 32.7 
days the error in recovering the initial positions of each of the planets is of 
order 10" 5 AU or about 1500km. Note that Jupiter has in this time travelled 
2.5 x 10 15 km. 

The error in the longitude of Jupiter can be estimated if we assume that 
the energy error is mainly in the the orbit of Jupiter. The relative energy 
error is proportional to the relative error in orbital frequency so the error in 
longitude is proportional to the integral of the relative energy error: AA « 
tnAE(t)/E, where n is the mean motion of Jupiter and t is the time of 
integration. Since the energy error grows linearly with time the position 
error grows with the square of the time. We find that the accumulated error 
in the longitude of Jupiter after 100 million years is only about 4 minutes of 
arc. This is to be compared with the 50 degree accumulated error estimated 
for our 200 million year integrations. The error in the longitude of Jupiter 
after the full 845 million years is about 5 degrees. 

We have directly measured the integration error in the determination of 
the position of Pluto by integrating forward and backward over intervals as 
long as 3 million years to determine how well we can reproduce the initial 
conditions. Over such short intervals the roundtrip error in the position of 
Pluto grows as a power of the time with an exponent near 2. The error in 
position is approximately 1.3 x 10" 19 * 2 AU (where t is in years). This growth 
of error is almost entirely in the integration of Pluto's orbit; the roundtrip 
error is roughly the same when we integrate the whole system and when we 
integrate Pluto in the field of the Sun only. It is interesting to note that 
in the integrations using the 32.7 day stepsize the position errors in all the 



planets are comparable. Extrapolation of the roundtrip error for Pluto over 
the full 845 million year integration gives an error in longitude of less than 
10 minutes of arc. 

Features of the orbital elements of Pluto 

Our new integrations confirm and extend our previous report of extremely 
long periods in the orbital elements of Pluto. 

For the first 450 million years of our integration we recorded the state 
of the system every 499983 days (about 1369 years) of simulated time. For 
the second 400 million years we sampled 16 times less frequently. For each 
planet the positions and velocities were converted to orbital elements relative 
to the center of mass of those bodies with smaller semimajor axes [9]. The 
elements were used to form the variables h, k, p, and q, defined as 

h = e sin w k = e cos w 

p = sin(e/2) sin ft q = sin(i/2) cos Q 

where i is the inclination, e is the eccentricity, Q, is the longitude of the 
ascending node, and w is the longitude of the perihelion. The variables p 
and q specify the orientation of the orbital plane, and the variables h and k 
specify the eccentricity and the orientation of the orbit in the orbital plane. 
h, k, p, and q are natural variables for looking at the long-term behavior of 
planetary systems. 

Figure 1 shows the variation of h with time. The largest component is 
the 3.7 million year regression of the longitude of perihelion. The 27 million 
year component we previously reported is clearly visible, as is the 137 million 
year component. The quantity p (not shown) has similar features. 

Figure 2 shows the inclination of Pluto as a function of time. Besides 
the major 3.8 million year component we can clearly discern the 34 million 
year component we previously reported. Although there is no continuing 
secular decline in the inclination, there is a component with a period near 150 
million years and evidence for a component with a period of approximately 
600 million years. The existence of significant orbital variations with such 
long periods is quite surprising. For quasiperiodic trajectories we expect 
to find frequencies which are low order combinations of a few fundamental 
frequencies (one per degree of freedom). The natural timescale for the long 
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Figure 1: The orbital element h = esinw for Pluto over 845 million years. 
On this scale the dominant period (the 3.7 million year circulation of the 
longitude of perihelion) is barely resolved. The most obvious component has 
a period of 137 million years. The sampling interval was increased in the 
second half of our integration. 



term evolution of a quasiperiodic planetary systems is set by the periods 
of the circulation of the nodes and perihelia. The presence of significant 
components with periods much longer than these indicates strong low-order 
resonances. An abundance of low-order resonances can give rise to chaotic 
behavior. Pluto is involved in at least five resonances, a sixth is possibly 
indicated by the 600 million year component that we can see in figure 2. Each 
time the motion of Pluto has been investigated over longer intervals longer 
periods have been found. This is suggestive that the motion is chaotic. 

Deterministic Chaotic Behavior 

In most conservative dynamical systems Newton's equations have both 
regular solutions and chaotic solutions. For some initial conditions the mo- 
tion is quasiperiodic; for others the motion is chaotic. Chaotic behavior is 
distinguished from quasiperiodic behavior by the way in which nearby trajec- 
tories diverge. [10,11] Nearby quasiperiodic trajectories diverge linearly with 
time on average, whereas nearby chaotic trajectories diverge exponentially 
with time. Quasiperiodic motion can be reduced to motion on a multidimen- 
sional torus; the frequency spectrum of quasiperiodic motion has as many 
independent frequencies as degrees of freedom. The frequency spectrum of 
chaotic motion is more complicated, usually appearing to have a broad-band 
component. 

The Lyapunov exponents measure the average rates of exponential diver- 
gence of nearby orbits. The Lyapunov exponents are limits for large time of 
the quantity 7 = ln(d/d )/(t—to), where d is the distance between the trajec- 
tory and an infinitesimally nearby test trajectory, and t is the time. For any 
particular trajectory of an n-dimensional system there can be n distinct Lya- 
punov exponents, depending on the phase space direction from the reference 
trajectory to the nearby trajectory. In Hamiltonian systems the Lyapunov 
exponents are paired; for each non-negative exponent there is a non-positive 
exponent with equal magnitude. Thus an m degree-of-freedom Hamiltonian 
system can have at most m positive exponents. For chaotic trajectories the 
largest Lyapunov exponent is positive; for quasiperiodic trajectories all of 
the Lyapunov exponents are zero. 

Lyapunov exponents can be estimated from the time evolution of the 
phase space distance between a reference trajectory and nearby test trajec- 
tories [10,12]. The most straightforward approach is to simply follow the 
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Figure 2: The inclination of Pluto over 845 million years. Besides the 34 
million year component and the 150 million year component there appears 
to be a component with a period near 600 million years. 



trajectories of a small cloud of particles started with nearby initial condi- 
tions. With a sufficiently long integration we can determine if the distances 
between the particles in the cloud diverge exponentially or linearly. If the 
divergence is exponential then for each pair of particles in the cloud we ob- 
tain an estimate of the largest Lyapunov exponent. With this method the 
trajectories eventually diverge so much that they no longer sample the same 
neighborhood of the phase space. We could fix this by periodically rescaling 
the cloud to be near the reference trajectory, but we can even more directly 
study the behavior of trajectories in the neighborhood of a reference trajec- 
tory by integrating the variational equations along with the reference tra- 
jectory. In particular, let y' = f (y) be an autonomous system of first-order 
ordinary differential equations and y(t) be the reference trajectory. We de- 
fine a phase-space variational trajectory y + Sy and note that 6y satisfies 
a linear system of first-order ordinary differential equations with coefficients 
that depend on y(t), 8y = J >6y where the elements of the Jacobian matrix 
are Jij = dfi/dyj. 

Lyapunov Exponent of Pluto 

We estimated the largest Lyapunov exponent of Pluto by both the vari- 
ational and the phase-space distance methods over the second half of our 
845 million year run. Figure 3 shows the divergence of the logarithm of the 
phase-space distance in a representative 2-particle experiment and the growth 
of the logarithm of the variational phase-space distance. We measured the 
phase-space distance by the ordinary Euclidean norm in the 6-dimensional 
space with position and velocity coordinates. We measured position in AU 
and velocity in AU/day. Since the magnitude of the velocity in these units is 
small compared to the magnitude of the position, the phase space distance is 
effectively equivalent to the positional distance, and we refer to phase space 
distances in terms of AU. For both traces in this plot the average growth 
is linear, indicating exponential divergence of nearby trajectories with an e- 
folding time of approximately 20 million years. The shapes of these graphs 
are remarkably similar until the two-particle divergence grows to about an 
astronomical unit, verifying that the motion in the neighborhood of Pluto is 
properly represented. A more conservative representation of this data is to 
plot the logarithm of 7 versus the logarithm of time, as shown in figure 4. 
The leveling off of this graph indicates a positive Lyapunov exponent. 
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Figure 3: The exponential divergence of nearby trajectories is indicated by 
the average linear growth of the logarithms of the distance measures as a 
function of time. In the upper trace we see the growth of the variational 
distance around a reference trajectory. In the lower trace we see how two 
Plutos diverge with time. The distance saturates near 80AU when the Plu- 
tos are on opposite sides of the Sun. The variational method of studying 
neighboring trajectories does not have the problem of saturation. Note that 
the two methods are in excellent agreement until the two-trajectory method 
has nearly saturated. 
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Figure 4: The conventional representation of the Lyapunov exponent calcu- 
lation, the logarithm of 7 versus the logarithm of time. Convergence to a 
positive exponent is indicated by a leveling off; for regular trajectories this 
plot approaches a line with slope minus one. 
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To study the details of the divergence of nearby trajectories we expand 
the early portion of the divergence graph. Figure 5 shows the logarithms 
of the distances between several pairs of particles in the cloud of Plutos in 
our integration versus the logarithm of time. The separation starts out as 
a power law with an exponent near 3/2. Only after some time does the 
exponential take off. The power law is dominated by the exponential only 
after the rate of growth of the exponential exceeds the rate of growth of 
the power law. This suggests that the portion of the divergence of nearby 
trajectories that results only from the numerical error fits a 3/2 power law 
and that this error "seeds" the exponential divergence that is the hallmark of 
chaos. We tested this hypothesis by integrating a cloud of test particles with 
the orbital elements of Pluto in the field of the Sun alone. The divergence of 
these Kepler "Plutos" grows as 3.16 x 10" 17 * 3 / 2 AU. This is identical to the 
initial divergence of the Plutos in the complete dynamical system, showing 
that two-body numerical error completely accounts for the initial divergence. 

Only the second half of the integration was used in the computation of 
the Lyapunov exponents, because the measurement in the first half of our 
integration was contaminated by over-vigorous application of the renormal- 
ization method, and gave a Lyapunov exponent about a factor of four too 
large. The renormalization interval was only 275,000 years, which was far too 
small. The renormalization interval must be long enough that the divergence 
of neighboring trajectories is dominated by the exponential divergence due 
to the sensitive dependence on initial conditions rather than the power law 
divergence due to the accumulation of numerical errors. In our experiment 
the renormalization interval should have been greater than 30 million years. 
It is important to emphasize that the variational method of measuring the 
Lyapunov exponent has none of these problems. 

Discussion 

Usually the measurement of a positive Lyapunov exponent provides a 
confirmation of what is already visible to the eye, chaotic trajectories usually 
look irregular. In this case the plots of Pluto's orbital elements do not look 
particularly irregular (see figures 1 and 2). The irregularity of the motion 
does manifest itself in the power spectra. For a quasiperiodic trajectory 
the power spectrum of any orbital element is composed of integral linear 
combinations of fundamental frequencies, where the number of fundamental 
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Figure 5: Common logarithm of the distance between several pairs of Plutos, 
in AU, versus the common logarithm of the time, in years. The initial seg- 
ment of the graph closely fits a 3/2 power law (dashed line). The solid line 
is an exponential chosen to fit the long-time divergence of Plutos. The expo- 
nential growth takes over when its slope exceeds the slope of the power-law. 
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frequencies is equal to the number of degrees of freedom. The power spectrum 
of chaotic trajectory is more complicated, usually appearing to have some 
broad-band component. 

A portion of the power spectrum of Neptune's h is shown in figure 6. For 
comparison the same portion of the power spectrum of Pluto's h is shown 
in figure 7. Hanning windows have been used to reduce spectral leakage; 
only the densely sampled part of the run was used in the computation of 
the Fourier transform. The spectrum of Neptune is quite complicated but 
there is no evidence that it is not a line spectrum. On the other hand the 
spectrum of Pluto does appear to have a broad-band component. Note that 
both of these spectra are computed from the same integration run, using the 
same numerical methods. They are subject to the same error processes, so 
the differences we see are dynamical in origin. The broad-band character of 
the Pluto spectrum is consistent with the chaotic nature of the motion as 
indicated by our measurement of a positive Lyapunov exponent. 

The lack of obvious irregularity in the orbital elements of Pluto indicates 
that the portion of the chaotic zone in which Pluto is currently moving is 
rather small. Since the global structure of the chaotic zone is not known it is 
not possible for us to predict whether more irregular motions are likely. If the 
small chaotic zone in which Pluto is found connects to a larger chaotic region 
relatively sudden transitions can be made to more irregular motion. This 
actually occurs for the motion of asteroids near the 3/1 Kirkwood gap [13]. 

On the other hand, the fact that the timescale for divergence is only an 
order of magnitude larger than the fundamental timescales of the system 
indicates that the chaotic behavior is robust. It is not a narrow chaotic zone 
associated with a high order resonance. Even though we do not know the 
sensitivity of the observed chaotic behavior to the uncertainties in parameters 
and initial conditions, and unmodelled effects the large Lyapunov exponent 
suggests that the chaotic behavior is characteristic of a range of solar systems 
including the actual solar system. 

Conclusion 

Our numerical model indicates that the motion of Pluto is chaotic. The 
largest Lyapunov exponent is about 10 _73 yr _1 . Thus the e- folding time for 
the divergence of trajectories is about 20 million years. This is a remarkably 
short e-folding time, considering the age of the Solar System. 
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Figure 6: A portion of the power spectrum of Neptune's h. The spectrum is 
quite complicated but there is no evidence of a broadband component; the 
spectrum is consistent with a line spectrum. 
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Figure 7: A portion of the power spectrum of Pluto's h. There appears 
to be a broad-band component to the spectrum. This is consistent with 
the chaotic character of the motion of Pluto as indicated by the positive 
Lyapunov exponent. 
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